Setup

Set up workspace, set options, and load required packages.

rm(list=ls(all=TRUE)) 
knitr::opts_chunk$set(root.dir = "~/R",warning=FALSE, message=FALSE)
library(gridExtra)
library(multcomp)
library(emmeans)
library(ggplot2)
library(dplyr)
library(effects)
library(plyr)
library(lattice)
library(glmmTMB)
library(effects)
library(lsmeans)
library(MuMIn)
library(MASS)
library(car)
library(FSA)
library(lmerTest)
library(lme4)
library(blmeco)
library(ggsci)
library(coxme)
library(survival)
library(cowplot)
library(RColorBrewer)
library(tidyverse)
library(partitionBEFsp)
library(factoextra)
library(ggfortify)
library(cowplot)
library(survminer)

1. GROW OUT SURVIVORSHIP

Load data

Survivorship is measured as binary on each individual spat in the experiment.

#load data
rearing<-read.csv("Data/GrowOut_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA")
rearing$Quadrant<-as.factor(rearing$Quadrant)
rearing$Diversity.Type<-relevel(rearing$Diversity.Type, ref="Individual")

#rearing$Total.Juveniles<-as.factor(rearing$Total.Juveniles)
#rearing$Fusion.Partners<-as.factor(rearing$Fusion.Partners)
#rearing$Richness<-as.factor(rearing$Richness)

Subset the final timepoint for additional analyses.

final<-rearing[which(rearing$Timepoint == "Final"),]
multiple<-rearing[which(rearing$Community == "Multiple"),]
multiple_final<-multiple[which(multiple$Timepoint == "Final"),]

Analysis

Analyze with a logistic binomial regression model.

Analyze at the final timepoint.

#with final dataset
survmod2<-glmer(Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial)
summary(survmod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +  
##     (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
##    Data: final
## 
##      AIC      BIC   logLik deviance df.resid 
##    224.9    255.5   -103.5    206.9      212 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.6175 -0.6293  0.3049  0.5167  1.4715 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 3.186e-01 0.5644532
##  Colony:Parent.Site (Intercept) 1.488e-02 0.1219734
##  Parent.Site        (Intercept) 1.634e-07 0.0004042
##  Tank               (Intercept) 1.709e-01 0.4133998
## Number of obs: 221, groups:  
## Slide, 115; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)   
## (Intercept)               1.450767   0.519055   2.795  0.00519 **
## Richness                 -0.007283   0.267024  -0.027  0.97824   
## CommunityMultiple        -1.847623   0.618356  -2.988  0.00281 **
## Fusion.Partners           3.036311   0.946317   3.209  0.00133 **
## Richness:Fusion.Partners -0.432505   0.319527  -1.354  0.17587   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rchnss CmmntM Fsn.Pr
## Richness    -0.503                     
## CmmntyMltpl -0.231 -0.556              
## Fusn.Prtnrs -0.105  0.414 -0.505       
## Rchnss:Fs.P  0.145 -0.474  0.405 -0.902
## convergence code: 0
## Model failed to converge with max|grad| = 0.00237872 (tol = 0.001, component 1)
Anova(survmod2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                            Chisq Df Pr(>Chisq)    
## Richness                  0.5759  1   0.447925    
## Community                 8.9279  1   0.002808 ** 
## Fusion.Partners          21.1261  1    4.3e-06 ***
## Richness:Fusion.Partners  1.8322  1   0.175871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#final dataset
summary(effect(c("Richness", "Fusion.Partners"), xlevels=4, survmod2))
## 
##  Richness*Fusion.Partners effect
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.5201784 0.9361096 0.9949754 0.9996265
##        2 0.5183603 0.9041990 0.9880627 0.9986242
##        3 0.5165418 0.8587551 0.9719082 0.9949464
##        4 0.5147229 0.7966036 0.9353234 0.9816172
## 
##  Lower 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.3287517 0.8160149 0.9478899 0.9855863
##        2 0.3643003 0.7954607 0.9392540 0.9827277
##        3 0.3149904 0.7068390 0.8811552 0.9520535
##        4 0.2323719 0.5215122 0.6410654 0.7004323
## 
##  Upper 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.7058580 0.9797581 0.9995363 0.9999905
##        2 0.6690068 0.9581692 0.9977482 0.9998920
##        3 0.7128515 0.9387678 0.9938440 0.9994880
##        4 0.7879786 0.9336587 0.9915323 0.9991807
summary(effect("Fusion.Partners", xlevels=4, survmod2))
## 
##  Fusion.Partners effect
## Fusion.Partners
##         0         1         2         3 
## 0.5186977 0.9110372 0.9898281 0.9989197 
## 
##  Lower 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.3642350 0.8027598 0.9426089 0.9840155 
## 
##  Upper 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.6696673 0.9626407 0.9982685 0.9999280
summary(effect("Community", xlevels=4, survmod2))
## 
##  Community effect
## Community
## Individual   Multiple 
##  0.9480743  0.7421160 
## 
##  Lower 95 Percent Confidence Limits
## Community
## Individual   Multiple 
##  0.8373723  0.5918728 
## 
##  Upper 95 Percent Confidence Limits
## Community
## Individual   Multiple 
##  0.9847894  0.8509763
dispersion_glmer(survmod2) #no overdispersion
## [1] 0.9357489
plot(residuals(survmod2)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)

Create plot for effect of community type.

eff.surv1 <- predictorEffect("Community", survmod2)
surv.effects1<-plot(eff.surv1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Number of Juveniles")), 
                   main=""); surv.effects1

                   #lattice=list(key.args=list(space="right",
                                              #border=FALSE, 
                                              #title=expression(bold("Genotypic \nRichness")),
                                              #cex=1, 
                                              #cex.title=1)));surv.effects1

Create plots for fusion partner effects.

survmod2a<-glmer(Survivorship ~ Richness + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))

eff.surv2 <- predictorEffect("Fusion.Partners", survmod2a, xlevels=4, rug=TRUE)
surv.effects2<-plot(eff.surv2,
                   lines=list(multiline=TRUE, 
                              col=c("lightblue","mediumblue", "blue4", "darkblue"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));surv.effects2

survmod2b<-glmer(Survivorship ~ Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))

eff.surv2b <- predictorEffect("Fusion.Partners", survmod2b, xlevels=4, rug=TRUE)
surv.effects2b<-plot(eff.surv2b,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="none",
                                              border=FALSE, 
                                              title=expression(bold("")),
                                              cex=1, 
                                              cex.title=1)));surv.effects2b

Generate mean values of survivorship for 1) number of juveniles and 2) within “multiple” juvenile groups, mean values for number of fusion partners.

meansurv <- ddply(final, c("Community"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv
##    Community   N      mean        sd         se max     lower    upper
## 1 Individual  58 0.7758621 0.4206554 0.05523477   1 0.3552066 1.196518
## 2   Multiple 163 0.6687117 0.4721270 0.03697984   1 0.1965847 1.140839
meansurv2 <- ddply(multiple_final, c("Fusion.Partners"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv2
##   Fusion.Partners  N      mean        sd         se max      lower    upper
## 1               0 75 0.3466667 0.4791133 0.05532324   1 -0.1324466 0.825780
## 2               1 48 0.9583333 0.2019409 0.02914766   1  0.7563924 1.160274
## 3               2 24 0.9166667 0.2823299 0.05763034   1  0.6343368 1.198997
## 4               3 16 0.9375000 0.2500000 0.06250000   1  0.6875000 1.187500

2. GROW OUT FUSION

Analysis

Analyze with a logistic binomial regression model.

Final timepoint dataset.

Analyze fusion as a product of genotypic richness within “multiple” juvenile groups.

#with final dataset
fusemod1<-glmer(Fusion ~ Richness + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))
summary(fusemod1)
Anova(fusemod1, type="II") # no effects on rate of fusion by the end point

dispersion_glmer(fusemod1) #no overdispersion
plot(residuals(fusemod1)) + abline(h=0, lty=2)

There is an effect of genotypic diversity on fusion rates. Plot the relationship.

eff.fuse1 <- predictorEffect("Richness", fusemod1, xlevels=4, rug=TRUE)
fuse.effects1<-plot(eff.fuse1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1),
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   xlab=expression(bold("Richness")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects1

Generate a summary of proportion fused at the end of the grow-out period.

meanfusion <- ddply(multiple_final, c("Richness"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanfusion
##   Richness  N      mean        sd         se max       lower     upper
## 1        1 59 0.5423729 0.5024778 0.06541704   1  0.03989509 1.0448507
## 2        2 52 0.7115385 0.4574670 0.06343925   1  0.25407150 1.1690054
## 3        3 28 0.1428571 0.3563483 0.06734350   1 -0.21349118 0.4992055
## 4        4 24 0.6250000 0.4945354 0.10094661   1  0.13046464 1.1195354

3. GROW OUT GROWTH

Analysis

Growth is measured as polyp expansion in each spat by number of spat in each sibling type and number of fused individuals. Growth is measured in final dataset as all spat started with 1 polyp each.

Analyze polyp expansion with the final dataset and present means.

polypmod1<-glmer(Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family = poisson)
summary(polypmod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +  
##     (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
##    Data: final
## 
##      AIC      BIC   logLik deviance df.resid 
##    590.6    618.1   -286.3    572.6      148 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1240 -0.5469 -0.1021  0.4563  1.4273 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 1.858e-09 4.311e-05
##  Colony:Parent.Site (Intercept) 7.237e-02 2.690e-01
##  Parent.Site        (Intercept) 1.445e-03 3.802e-02
##  Tank               (Intercept) 2.789e-03 5.282e-02
## Number of obs: 157, groups:  
## Slide, 92; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.279359   0.147887   8.651   <2e-16 ***
## Richness                 -0.021274   0.091952  -0.231    0.817    
## CommunityMultiple        -0.095335   0.160455  -0.594    0.552    
## Fusion.Partners           0.068683   0.143007   0.480    0.631    
## Richness:Fusion.Partners -0.009974   0.052539  -0.190    0.849    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rchnss CmmntM Fsn.Pr
## Richness    -0.623                     
## CmmntyMltpl  0.044 -0.611              
## Fusn.Prtnrs -0.291  0.579 -0.673       
## Rchnss:Fs.P  0.384 -0.723  0.582 -0.912
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(polypmod1, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                           Chisq Df Pr(>Chisq)
## Richness                 0.2843  1     0.5939
## Community                0.3530  1     0.5524
## Fusion.Partners          0.5587  1     0.4548
## Richness:Fusion.Partners 0.0360  1     0.8494
qqPlot(residuals(polypmod1))

## 389 393 
##  73  74
meangrowth <- ddply(final, c("Community"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meangrowth
##    Community   N     mean       sd        se max    lower    upper
## 1 Individual  45 3.533333 1.778661 0.2651472   8 1.754672 5.311995
## 2   Multiple 112 3.357143 1.558800 0.1472928   7 1.798343 4.915943

Plot the model outputs.

eff.polyps1 <- predictorEffect("Fusion.Partners", polypmod1, xlevels=4, rug=TRUE)
polyp.effects1<-plot(eff.polyps1,
                   lines=list(multiline=TRUE, 
                              col=c("lightblue","mediumblue", "blue4", "darkblue"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("Total Polyp Expansion")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));polyp.effects1

Combine figures for Grow-Out Period.

grow_out_figs<-plot_grid(surv.effects1, surv.effects2b, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(0.8,1), label_size = 20, label_y=1, align="h");grow_out_figs

ggsave(filename="Figures/grow_out_figs.png", plot=grow_out_figs, dpi=500, width=10, height=6, units="in")

4. STRESS SURVIVORSHIP

Load data

Survivorship is measured as binary on each individual spat in the experiment.

#load data
stress<-read.csv("Data/Stress_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA")
stress$Diversity.Type<-relevel(stress$Diversity.Type, ref="Individual")
stress$Polyps<-as.numeric(stress$Polyps)

Subset the final timepoint, “multiple” juvenile datasets for additional analyses.

final_stress<-stress[which(stress$Timepoint == "End"),]
multiple_stress<-stress[which(stress$Community == "Multiple"),]
multiple_final_stress<-multiple_stress[which(multiple_stress$Timepoint == "End"),]
d26_stress<-stress[which(stress$Timepoint == "S2D11"),]
d23_stress<-stress[which(stress$Timepoint == "S2D8"),]

Analysis

Analyze with a logistic binomial regression model.

Check for effect of community type (“individual” or “multiple” juveniles) on survivorship and display means.

#with full time dataset
survmod3b<-glmer(Survivorship ~ Day * Treatment * Community + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit) 

summary(survmod3b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Survivorship ~ Day * Treatment * Community + (1 | Parent.Site/Colony) +  
##     (1 | Treatment:Tank) + (1 | Slide)
##    Data: stress
## 
##      AIC      BIC   logLik deviance df.resid 
##    449.0    508.9   -212.5    425.0     1076 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -110.079    0.001    0.028    0.118    5.533 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 2.605e+00 1.614e+00
##  Colony:Parent.Site (Intercept) 5.717e-01 7.561e-01
##  Treatment:Tank     (Intercept) 5.061e-09 7.114e-05
##  Parent.Site        (Intercept) 1.645e-08 1.282e-04
## Number of obs: 1088, groups:  
## Slide, 90; Colony:Parent.Site, 20; Treatment:Tank, 9; Parent.Site, 6
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          6.7660321  1.7130404   3.950 7.82e-05 ***
## Day                                 -0.0835370  0.0592967  -1.409 0.158896    
## TreatmentHigh                        0.8459104  1.9259267   0.439 0.660500    
## CommunityMultiple                    3.4775795  2.4771976   1.404 0.160368    
## Day:TreatmentHigh                   -0.2886216  0.0823811  -3.503 0.000459 ***
## Day:CommunityMultiple               -0.1855821  0.0878809  -2.112 0.034708 *  
## TreatmentHigh:CommunityMultiple      2.8355215  3.2070137   0.884 0.376608    
## Day:TreatmentHigh:CommunityMultiple  0.0001848  0.1226640   0.002 0.998798    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Day    TrtmnH CmmntM Dy:TrH Dy:CmM TrH:CM
## Day         -0.835                                          
## TreatmntHgh -0.700  0.666                                   
## CmmntyMltpl -0.586  0.538  0.484                            
## Dy:TrtmntHg  0.397 -0.633 -0.837 -0.348                     
## Dy:CmmntyMl  0.507 -0.652 -0.450 -0.917  0.471              
## TrtmntHg:CM  0.412 -0.395 -0.600 -0.760  0.504  0.702       
## Dy:TrtmH:CM -0.295  0.436  0.563  0.641 -0.648 -0.713 -0.924
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(survmod3b, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                           Chisq Df Pr(>Chisq)    
## Day                     65.5353  1  5.708e-16 ***
## Treatment               37.2722  1  1.027e-09 ***
## Community                3.3205  1   0.068421 .  
## Day:Treatment           21.1660  1  4.212e-06 ***
## Day:Community            9.0584  1   0.002615 ** 
## Treatment:Community      5.3796  1   0.020374 *  
## Day:Treatment:Community  0.0000  1   0.998798    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3b) #no overdispersion
## [1] 0.5599041
plot(residuals(survmod3b)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)
meanstresssurv1 <- ddply(final_stress, c("Community", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv1
##    Community Treatment  N       mean        sd         se max      lower
## 1 Individual   Ambient 21 0.95238095 0.2182179 0.04761905   1  0.7341631
## 2 Individual      High 24 0.00000000 0.0000000 0.00000000   0  0.0000000
## 3   Multiple   Ambient 46 0.76086957 0.4312660 0.06358670   1  0.3296036
## 4   Multiple      High 65 0.04615385 0.2114510 0.02622727   1 -0.1652972
##       upper
## 1 1.1705988
## 2 0.0000000
## 3 1.1921355
## 4 0.2576049

Analyze within “multiple” groups as survival was significantly different between single and multiple juvenile slides.

#with full time dataset

survmod3<-glmer(Survivorship ~ Day + Treatment + Fusion.Partners + Richness + Day:Treatment + Day:Fusion.Partners + Day:Richness + Day:Treatment:Fusion.Partners + Day:Treatment:Richness + Day:Fusion.Partners:Richness:Treatment + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, subset=c(Community=="Multiple")) 

#fixed_form<-glm(Survivorship ~ Day * Treatment * Fusion.Partners *Richness, data=stress, family=binomial, subset=c(Community=="Multiple"))

#library("brglm2"); glm(fixed_form, data = stress, family = binomial, method="detect_separation")

summary(survmod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Survivorship ~ Day + Treatment + Fusion.Partners + Richness +  
##     Day:Treatment + Day:Fusion.Partners + Day:Richness + Day:Treatment:Fusion.Partners +  
##     Day:Treatment:Richness + Day:Fusion.Partners:Richness:Treatment +  
##     (1 | Parent.Site/Colony) + (1 | Treatment:Tank) + (1 | Slide/Sample)
##    Data: stress
##  Subset: c(Community == "Multiple")
## 
##      AIC      BIC   logLik deviance df.resid 
##    222.5    301.5    -94.2    188.5      756 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.6798  0.0000  0.0000  0.0000  0.5145 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample:Slide       (Intercept) 940.14567 30.6618 
##  Slide              (Intercept)   0.04049  0.2012 
##  Colony:Parent.Site (Intercept)   0.01752  0.1324 
##  Treatment:Tank     (Intercept)   8.33415  2.8869 
##  Parent.Site        (Intercept)   0.03312  0.1820 
## Number of obs: 773, groups:  
## Sample:Slide, 111; Slide, 45; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                                Estimate Std. Error z value
## (Intercept)                                    -2.29443   39.97253  -0.057
## Day                                             0.55107    1.24309   0.443
## TreatmentHigh                                 111.46301   35.69621   3.123
## Fusion.Partners                               100.44513   51.46355   1.952
## Richness                                       22.39796   34.13494   0.656
## Day:TreatmentHigh                              -4.98049    1.87839  -2.651
## Day:Fusion.Partners                            -2.73031    1.88657  -1.447
## Day:Richness                                   -0.90135    1.21483  -0.742
## Day:TreatmentHigh:Fusion.Partners              -1.42109    0.99729  -1.425
## Day:TreatmentHigh:Richness                     -0.02063    0.80677  -0.026
## Day:TreatmentAmbient:Fusion.Partners:Richness  -0.09927    0.27596  -0.360
## Day:TreatmentHigh:Fusion.Partners:Richness      0.07110    0.22617   0.314
##                                               Pr(>|z|)   
## (Intercept)                                    0.95423   
## Day                                            0.65754   
## TreatmentHigh                                  0.00179 **
## Fusion.Partners                                0.05097 . 
## Richness                                       0.51172   
## Day:TreatmentHigh                              0.00801 **
## Day:Fusion.Partners                            0.14783   
## Day:Richness                                   0.45811   
## Day:TreatmentHigh:Fusion.Partners              0.15417   
## Day:TreatmentHigh:Richness                     0.97959   
## Day:TreatmentAmbient:Fusion.Partners:Richness  0.71905   
## Day:TreatmentHigh:Fusion.Partners:Richness     0.75326   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Day    TrtmnH Fsn.Pr Rchnss Dy:TrH Dy:F.P Dy:Rch Dy:TH:F.P
## Day         -0.633                                                           
## TreatmntHgh  0.030 -0.317                                                    
## Fusn.Prtnrs  0.516 -0.509  0.005                                             
## Richness    -0.917  0.826 -0.219 -0.590                                      
## Dy:TrtmntHg -0.392  0.059 -0.623 -0.043  0.331                               
## Dy:Fsn.Prtn -0.611  0.393  0.087 -0.936  0.590  0.133                        
## Day:Richnss  0.627 -0.990  0.337  0.519 -0.838 -0.074 -0.404                 
## Dy:TrtH:F.P  0.061  0.336 -0.155 -0.324  0.137 -0.338  0.071 -0.337          
## Dy:TrtmnH:R  0.478  0.213 -0.071  0.098 -0.287 -0.680 -0.275 -0.210  0.453   
## Dy:TA:F.P:R  0.424  0.315 -0.245  0.088 -0.181 -0.371 -0.391 -0.313  0.552   
## Dy:TH:F.P:R  0.162 -0.148 -0.066  0.147 -0.178  0.307 -0.142  0.151 -0.652   
##             D:TH:R D:TA:F
## Day                      
## TreatmntHgh              
## Fusn.Prtnrs              
## Richness                 
## Dy:TrtmntHg              
## Dy:Fsn.Prtn              
## Day:Richnss              
## Dy:TrtH:F.P              
## Dy:TrtmnH:R              
## Dy:TA:F.P:R  0.752       
## Dy:TH:F.P:R -0.255  0.031
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 6 negative eigenvalues
Anova(survmod3, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                                          Chisq Df Pr(>Chisq)    
## Day                                     8.3067  1    0.00395 ** 
## Treatment                               5.5318  1    0.01867 *  
## Fusion.Partners                         0.6676  1    0.41390    
## Richness                                0.2974  1    0.58553    
## Day:Treatment                          16.8143  1  4.122e-05 ***
## Day:Fusion.Partners                     0.8008  1    0.37084    
## Day:Richness                            1.0944  1    0.29549    
## Day:Treatment:Fusion.Partners           4.1122  1    0.04257 *  
## Day:Treatment:Richness                  0.3144  1    0.57499    
## Day:Treatment:Fusion.Partners:Richness  0.2355  2    0.88894    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3) #no overdispersion
## [1] 0.2731621
plot(residuals(survmod3)) + abline(h=0, lty=2) #Dispersed randomly

## integer(0)
meanstresssurv2 <- ddply(final_stress, c("Fusion.Partners", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv2
##   Fusion.Partners Treatment  N       mean        sd         se max      lower
## 1               0   Ambient 30 0.83333333 0.3790490 0.06920457   1  0.4542843
## 2               0      High 50 0.04000000 0.1979487 0.02799417   1 -0.1579487
## 3               1   Ambient 16 0.81250000 0.4031129 0.10077822   1  0.4093871
## 4               1      High 15 0.06666667 0.2581989 0.06666667   1 -0.1915322
## 5               2   Ambient  9 0.77777778 0.4409586 0.14698618   1  0.3368192
## 6               2      High 16 0.00000000 0.0000000 0.00000000   0  0.0000000
## 7               3   Ambient 12 0.83333333 0.3892495 0.11236664   1  0.4440839
## 8               3      High  8 0.00000000 0.0000000 0.00000000   0  0.0000000
##       upper
## 1 1.2123824
## 2 0.2379487
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000

Plot effect of number of juveniles over time.

stress$group<-paste(stress$Treatment, "-", stress$Community)

survmod3b2<-glmer(Survivorship ~ group * Day + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit) 

eff.surv3b2 <- predictorEffect(c("Day"), survmod3b2)
surv.effects3b2<-plot(eff.surv3b2,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(1, 2, 1, 2)), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Number of Juveniles")),
                                              cex=1, 
                                              cex.title=1)));surv.effects3b2

Examine curve estimates to identify 50% survivorship.

summary(effect(c("group", "Day"), xlevels=50, survmod3b2))
## 
##  group*Day effect
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9989505 0.9988891 0.9988237 0.9987553 0.9986828
##   Ambient - Multiple   0.9999639 0.9999570 0.9999487 0.9999389 0.9999273
##   High - Individual    0.9995124 0.9993782 0.9992059 0.9989885 0.9987118
##   High - Multiple      0.9999991 0.9999986 0.9999980 0.9999972 0.9999960
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9986050 0.9985238 0.9984379 0.9983470 0.9982494
##   Ambient - Multiple   0.9999132 0.9998966 0.9998769 0.9998534 0.9998250
##   High - Individual    0.9983533 0.9979030 0.9973299 0.9966008 0.9956574
##   High - Multiple      0.9999942 0.9999916 0.9999880 0.9999828 0.9999752
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9981475 0.9980398 0.9979240 0.9978033 0.9976756
##   Ambient - Multiple   0.9997916 0.9997518 0.9997037 0.9996472 0.9995799
##   High - Individual    0.9944741 0.9929706 0.9910287 0.9885985 0.9855198
##   High - Multiple      0.9999644 0.9999489 0.9999263 0.9998942 0.9998482
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9975384 0.9974066 0.9972439 0.9970711 0.9969144
##   Ambient - Multiple   0.9994984 0.9994107 0.9992889 0.9991419 0.9989919
##   High - Individual    0.9815577 0.9770442 0.9704068 0.9619250 0.9528374
##   High - Multiple      0.9997810 0.9996944 0.9995493 0.9993353 0.9990726
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9967210 0.9965456 0.9963292 0.9961329 0.9958907
##   Ambient - Multiple   0.9987836 0.9985711 0.9982760 0.9979751 0.9975572
##   High - Individual    0.9396333 0.9256365 0.9055717 0.8846468 0.8552516
##   High - Multiple      0.9986325 0.9980927 0.9971889 0.9960815 0.9942302
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9956711 0.9954001 0.9951544 0.9948512 0.9945764
##   Ambient - Multiple   0.9971311 0.9965396 0.9959368 0.9951002 0.9942481
##   High - Individual    0.8253269 0.7844978 0.7443186 0.6916290 0.6420357
##   High - Multiple      0.9919665 0.9881940 0.9836013 0.9759955 0.9668167
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9942373 0.9939299 0.9935506 0.9931477 0.9927826
##   Ambient - Multiple   0.9930663 0.9918634 0.9901965 0.9881921 0.9861554
##   High - Individual    0.5801574 0.5249503 0.4598598 0.3961105 0.3440645
##   High - Multiple      0.9518077 0.9340054 0.9056046 0.8667256 0.8233276
##                       Day
## group                       22.9      23.5      24.2      24.8      25.5
##   Ambient - Individual 0.9923322 0.9919240 0.9914204 0.9909641 0.9904013
##   Ambient - Multiple   0.9833386 0.9804813 0.9765376 0.9725465 0.9670535
##   High - Individual    0.2878140 0.2442424 0.1993517 0.1660498 0.1330014
##   High - Multiple      0.7595582 0.6935999 0.6054450 0.5237196 0.4270620
##                       Day
## group                       26.1       26.8       27.4       28.1      28.7
##   Ambient - Individual 0.9898913 0.98926241 0.98869264 0.98799004 0.9873536
##   Ambient - Multiple   0.9615127 0.95391658 0.94628912 0.93588877 0.9255105
##   High - Individual    0.1092707 0.08635276 0.07027067 0.05502712 0.0444949
##   High - Multiple      0.3481675 0.26582641 0.20600815 0.14957280 0.1119266
##                       Day
## group                        29.4         30       30.7       31.3         32
##   Ambient - Individual 0.98656895 0.98585830 0.98498226 0.98418900 0.98321132
##   Ambient - Multiple   0.91146307 0.89756352 0.87893523 0.86070948 0.83659930
##   High - Individual    0.03463450 0.02789024 0.02162627 0.01736948 0.01343575
##   High - Multiple      0.07870967 0.05768922 0.03984634 0.02887951 0.01976041
## 
##  Lower 95 Percent Confidence Limits
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9688623 0.9690922 0.9693058 0.9694989 0.9696732
##   Ambient - Multiple   0.9981187 0.9979334 0.9977283 0.9975049 0.9972591
##   High - Individual    0.9920845 0.9906655 0.9889803 0.9870133 0.9846947
##   High - Multiple      0.9999701 0.9999606 0.9999479 0.9999313 0.9999094
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9698303 0.9699644 0.9700772 0.9701676 0.9702353
##   Ambient - Multiple   0.9969842 0.9966859 0.9963574 0.9959955 0.9955902
##   High - Individual    0.9819168 0.9786896 0.9748889 0.9704148 0.9650636
##   High - Multiple      0.9998801 0.9998418 0.9997913 0.9997246 0.9996350
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9702769 0.9702921 0.9702789 0.9702358 0.9701607
##   Ambient - Multiple   0.9951498 0.9946641 0.9941193 0.9935266 0.9928718
##   High - Individual    0.9588613 0.9515783 0.9428952 0.9328710 0.9211579
##   High - Multiple      0.9995181 0.9993637 0.9991558 0.9988845 0.9985255
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9700495 0.9699156 0.9697168 0.9694682 0.9692122
##   Ambient - Multiple   0.9921364 0.9913992 0.9904473 0.9893846 0.9883749
##   High - Individual    0.9072779 0.8926825 0.8729571 0.8499648 0.8273689
##   High - Multiple      0.9980419 0.9974650 0.9965725 0.9953636 0.9939910
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9688590 0.9685055 0.9680282 0.9675585 0.9669328
##   Ambient - Multiple   0.9870676 0.9858233 0.9842088 0.9826685 0.9806649
##   High - Individual    0.7973226 0.7681915 0.7300762 0.6938098 0.6473968
##   High - Multiple      0.9918652 0.9894505 0.9857101 0.9814630 0.9748907
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9663236 0.9655190 0.9647409 0.9637193 0.9627359
##   Ambient - Multiple   0.9787481 0.9762472 0.9738466 0.9707029 0.9676733
##   High - Individual    0.6043419 0.5508357 0.5028190 0.4453401 0.3958471
##   High - Multiple      0.9674434 0.9559596 0.9430171 0.9232269 0.9011890
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9614496 0.9602154 0.9586051 0.9567911 0.9550561
##   Ambient - Multiple   0.9636884 0.9598297 0.9547280 0.9488733 0.9431543
##   High - Individual    0.3392213 0.2927669 0.2422482 0.1967494 0.1622952
##   High - Multiple      0.8680761 0.8320749 0.7797593 0.7153722 0.6506739
##                       Day
## group                       22.9      23.5       24.2       24.8      25.5
##   Ambient - Individual 0.9527987 0.9506428 0.94784145 0.94516956 0.9417021
##   Ambient - Multiple   0.9355228 0.9280213 0.91794801 0.90798712 0.8945388
##   High - Individual    0.1276564 0.1026581 0.07860727 0.06193355 0.0464441
##   High - Multiple      0.5657819 0.4877617 0.39554313 0.31998675 0.2407995
##                       Day
## group                        26.1       26.8       27.4       28.1       28.7
##   Ambient - Individual 0.93839918 0.93411856 0.93004678 0.92477792 0.91977455
##   Ambient - Multiple   0.88118186 0.86309627 0.84511815 0.82082232 0.79679637
##   High - Individual    0.03603222 0.02660821 0.02041373 0.01490956 0.01134799
##   High - Multiple      0.18322882 0.12923146 0.09371102 0.06306962 0.04427288
##                       Day
## group                         29.4          30        30.7        31.3
##   Ambient - Individual 0.913312333 0.907188427 0.899297391 0.891838810
##   Ambient - Multiple   0.764632262 0.733264999 0.692078154 0.652889335
##   High - Individual    0.008223676 0.006223925 0.004485222 0.003380753
##   High - Multiple      0.028915333 0.019893620 0.012757874 0.008671295
##                       Day
## group                           32
##   Ambient - Individual 0.882256056
##   Ambient - Multiple   0.602984743
##   High - Individual    0.002426434
##   High - Multiple      0.005499545
## 
##  Upper 95 Percent Confidence Limits
##                       Day
## group                          0     0.653      1.31      1.96      2.61
##   Ambient - Individual 0.9999657 0.9999612 0.9999562 0.9999506 0.9999444
##   Ambient - Multiple   0.9999993 0.9999991 0.9999988 0.9999985 0.9999981
##   High - Individual    0.9999702 0.9999589 0.9999433 0.9999221 0.9998930
##   High - Multiple      1.0000000 1.0000000 0.9999999 0.9999999 0.9999998
##                       Day
## group                       3.27      3.92      4.57      5.22      5.88
##   Ambient - Individual 0.9999373 0.9999294 0.9999207 0.9999109 0.9998998
##   Ambient - Multiple   0.9999975 0.9999968 0.9999959 0.9999947 0.9999931
##   High - Individual    0.9998523 0.9997972 0.9997218 0.9996186 0.9994748
##   High - Multiple      0.9999997 0.9999996 0.9999993 0.9999989 0.9999983
##                       Day
## group                       6.53      7.18      7.84      8.49      9.14
##   Ambient - Individual 0.9998876 0.9998740 0.9998587 0.9998420 0.9998236
##   Ambient - Multiple   0.9999911 0.9999885 0.9999851 0.9999809 0.9999754
##   High - Individual    0.9992809 0.9990161 0.9986487 0.9981550 0.9974840
##   High - Multiple      0.9999974 0.9999959 0.9999936 0.9999900 0.9999844
##                       Day
## group                        9.8      10.4      11.1      11.8      12.4
##   Ambient - Individual 0.9998028 0.9997821 0.9997555 0.9997261 0.9996985
##   Ambient - Multiple   0.9999682 0.9999599 0.9999475 0.9999313 0.9999134
##   High - Individual    0.9965576 0.9954292 0.9936503 0.9912023 0.9883944
##   High - Multiple      0.9999755 0.9999632 0.9999409 0.9999050 0.9998575
##                       Day
## group                       13.1      13.7      14.4        15      15.7
##   Ambient - Individual 0.9996634 0.9996306 0.9995892 0.9995507 0.9995024
##   Ambient - Multiple   0.9998868 0.9998576 0.9998142 0.9997666 0.9996959
##   High - Individual    0.9840224 0.9790595 0.9714308 0.9629020 0.9500353
##   High - Multiple      0.9997714 0.9996576 0.9994521 0.9991813 0.9986941
##                       Day
## group                       16.3        17      17.6      18.3      18.9
##   Ambient - Individual 0.9994579 0.9994024 0.9993517 0.9992890 0.9992323
##   Ambient - Multiple   0.9996189 0.9995047 0.9993806 0.9991973 0.9989992
##   High - Individual    0.9359643 0.9152970 0.8933854 0.8623573 0.8307878
##   High - Multiple      0.9980548 0.9969114 0.9954211 0.9927783 0.9893703
##                       Day
## group                       19.6      20.2      20.9      21.6      22.2
##   Ambient - Individual 0.9991628 0.9991006 0.9990252 0.9989470 0.9988782
##   Ambient - Multiple   0.9987079 0.9983946 0.9979371 0.9973571 0.9967406
##   High - Individual    0.7881158 0.7468255 0.6939347 0.6372250 0.5868072
##   High - Multiple      0.9834108 0.9758592 0.9629580 0.9439063 0.9210072
##                       Day
## group                       22.9      23.5      24.2      24.8      25.5
##   Ambient - Individual 0.9987962 0.9987249 0.9986410 0.9985688 0.9984850
##   Ambient - Multiple   0.9958518 0.9949165 0.9935835 0.9921980 0.9902508
##   High - Individual    0.5274205 0.4772442 0.4208513 0.3751904 0.3257642
##   High - Multiple      0.8845106 0.8432974 0.7825316 0.7198529 0.6365925
##                       Day
## group                       26.1      26.8      27.4      28.1      28.7
##   Ambient - Individual 0.9984139 0.9983324 0.9982640 0.9981867 0.9981227
##   Ambient - Multiple   0.9882570 0.9855000 0.9827247 0.9789551 0.9752286
##   High - Individual    0.2870447 0.2463006 0.2151503 0.1830334 0.1589008
##   High - Multiple      0.5598160 0.4690315 0.3943239 0.3148500 0.2553418
##                       Day
## group                       29.4        30       30.7       31.3         32
##   Ambient - Individual 0.9980511 0.9979928 0.99792836 0.99787649 0.99782004
##   Ambient - Multiple   0.9702584 0.9654314 0.95910208 0.95305189 0.94523397
##   High - Individual    0.1343732 0.1161636 0.09783686 0.08434176 0.07084932
##   High - Multiple      0.1968694 0.1558722 0.11759973 0.09182032 0.06845565

Create plot for effect of temperature on survival in multiple juvenile groups across number of fusion partners.

stress$group<-paste(stress$Treatment, "-", stress$Fusion.Partners)

survmod3a<-glmer(Survivorship ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial)

summary(effect(c("Day", "group"), xlevels=50, survmod3a))
## 
##  Day*group effect
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3   High - 0     High - 1
##   0       0.9996515   0.9999990   0.9999988   0.9999991 0.99974119 0.9999999995
##   0.653   0.9996149   0.9999987   0.9999985   0.9999989 0.99967271 0.9999999990
##   1.31    0.9995742   0.9999983   0.9999981   0.9999986 0.99958552 0.9999999982
##   1.96    0.9995297   0.9999979   0.9999975   0.9999983 0.99947644 0.9999999967
##   2.61    0.9994806   0.9999973   0.9999969   0.9999978 0.99933866 0.9999999939
##   3.27    0.9994255   0.9999965   0.9999960   0.9999973 0.99916166 0.9999999888
##   3.92    0.9993654   0.9999956   0.9999949   0.9999966 0.99894113 0.9999999796
##   4.57    0.9992992   0.9999944   0.9999935   0.9999957 0.99866268 0.9999999629
##   5.22    0.9992260   0.9999929   0.9999917   0.9999945 0.99831112 0.9999999324
##   5.88    0.9991438   0.9999909   0.9999895   0.9999931 0.99785968 0.9999998757
##   6.53    0.9990544   0.9999884   0.9999866   0.9999912 0.99729760 0.9999997737
##   7.18    0.9989557   0.9999853   0.9999830   0.9999890 0.99658842 0.9999995879
##   7.84    0.9988449   0.9999812   0.9999782   0.9999860 0.99567848 0.9999992427
##   8.49    0.9987243   0.9999762   0.9999723   0.9999824 0.99454672 0.9999986208
##   9.14    0.9985911   0.9999697   0.9999648   0.9999778 0.99312061 0.9999974885
##   9.8     0.9984417   0.9999614   0.9999551   0.9999718 0.99129383 0.9999953841
##   10.4    0.9982922   0.9999518   0.9999439   0.9999651 0.98921987 0.9999919730
##   11.1    0.9980996   0.9999377   0.9999273   0.9999552 0.98617675 0.9999846928
##   11.8    0.9978853   0.9999193   0.9999058   0.9999425 0.98228997 0.9999708098
##   12.4    0.9976825   0.9998994   0.9998824   0.9999288 0.97811848 0.9999492398
##   13.1    0.9974213   0.9998697   0.9998477   0.9999085 0.97203016 0.9999032059
##   13.7    0.9971741   0.9998375   0.9998098   0.9998867 0.96552684 0.9998316884
##   14.4    0.9968558   0.9997897   0.9997536   0.9998545 0.95609229 0.9996790825
##   15      0.9965547   0.9997377   0.9996924   0.9998197 0.94608877 0.9994420611
##   15.7    0.9961668   0.9996605   0.9996015   0.9997686 0.93171070 0.9989365612
##   16.3    0.9957999   0.9995766   0.9995024   0.9997133 0.91663663 0.9981521464
##   17      0.9953275   0.9994521   0.9993554   0.9996319 0.89527321 0.9964820804
##   17.6    0.9948807   0.9993166   0.9991953   0.9995440 0.87325018 0.9938982477
##   18.3    0.9943055   0.9991157   0.9989577   0.9994146 0.84267645 0.9884281352
##   18.9    0.9937616   0.9988971   0.9986989   0.9992748 0.81191650 0.9800473070
##   19.6    0.9930615   0.9985730   0.9983149   0.9990690 0.77043751 0.9626269868
##   20.2    0.9923997   0.9982206   0.9978969   0.9988469 0.73007757 0.9367555273
##   20.9    0.9915481   0.9976982   0.9972768   0.9985199 0.67771441 0.8859370088
##   21.6    0.9906020   0.9970230   0.9964745   0.9981004 0.62047378 0.8028772805
##   22.2    0.9897080   0.9962891   0.9956020   0.9976476 0.56851442 0.7007921559
##   22.9    0.9885583   0.9952024   0.9943092   0.9969815 0.50601579 0.5512081429
##   23.5    0.9874724   0.9940225   0.9929045   0.9962630 0.45222054 0.4139283323
##   24.2    0.9860764   0.9922772   0.9908260   0.9952069 0.39092373 0.2702669052
##   24.8    0.9847586   0.9903848   0.9885716   0.9940685 0.34091977 0.1755820234
##   25.5    0.9830654   0.9875907   0.9852427   0.9923969 0.28681031 0.1004629759
##   26.1    0.9814679   0.9845677   0.9816416   0.9905973 0.24477183 0.0603473376
##   26.8    0.9794168   0.9801171   0.9763424   0.9879594 0.20126260 0.0325806638
##   27.4    0.9774830   0.9753189   0.9706342   0.9851252 0.16879538 0.0189984697
##   28.1    0.9750021   0.9682869   0.9622797   0.9809815 0.13635299 0.0100534491
##   28.7    0.9726651   0.9607483   0.9533406   0.9765440 0.11287708 0.0058060161
##   29.4    0.9696698   0.9497787   0.9403684   0.9700824 0.09001823 0.0030530475
##   30      0.9668513   0.9381216   0.9266326   0.9631971 0.07383774 0.0017579282
##   30.7    0.9632430   0.9213463   0.9069597   0.9532348 0.05836465 0.0009226114
##   31.3    0.9598521   0.9037583   0.8864575   0.9427013 0.04757632 0.0005307535
##   32      0.9555170   0.8788699   0.8576668   0.9276079 0.03738431 0.0002783921
##        group
## Day         High - 2     High - 3
##   0     1.000000e+00 1.000000e+00
##   0.653 1.000000e+00 1.000000e+00
##   1.31  1.000000e+00 1.000000e+00
##   1.96  1.000000e+00 1.000000e+00
##   2.61  1.000000e+00 1.000000e+00
##   3.27  1.000000e+00 1.000000e+00
##   3.92  1.000000e+00 1.000000e+00
##   4.57  1.000000e+00 1.000000e+00
##   5.22  1.000000e+00 1.000000e+00
##   5.88  1.000000e+00 1.000000e+00
##   6.53  1.000000e+00 1.000000e+00
##   7.18  1.000000e+00 1.000000e+00
##   7.84  1.000000e+00 1.000000e+00
##   8.49  1.000000e+00 1.000000e+00
##   9.14  1.000000e+00 1.000000e+00
##   9.8   1.000000e+00 1.000000e+00
##   10.4  1.000000e+00 1.000000e+00
##   11.1  1.000000e+00 1.000000e+00
##   11.8  1.000000e+00 1.000000e+00
##   12.4  1.000000e+00 1.000000e+00
##   13.1  1.000000e+00 1.000000e+00
##   13.7  1.000000e+00 1.000000e+00
##   14.4  1.000000e+00 1.000000e+00
##   15    1.000000e+00 9.999999e-01
##   15.7  1.000000e+00 9.999997e-01
##   16.3  9.999999e-01 9.999993e-01
##   17    9.999996e-01 9.999978e-01
##   17.6  9.999988e-01 9.999941e-01
##   18.3  9.999956e-01 9.999816e-01
##   18.9  9.999862e-01 9.999515e-01
##   19.6  9.999476e-01 9.998494e-01
##   20.2  9.998352e-01 9.996025e-01
##   20.9  9.993730e-01 9.987669e-01
##   21.6  9.976180e-01 9.961820e-01
##   22.2  9.925485e-01 9.899812e-01
##   22.9  9.722222e-01 9.695402e-01
##   23.5  9.175687e-01 9.233977e-01
##   24.2  7.452137e-01 7.952102e-01
##   24.8  4.819230e-01 5.952330e-01
##   25.5  1.964146e-01 3.214383e-01
##   26.1  7.212871e-02 1.521096e-01
##   26.8  2.001693e-02 5.463164e-02
##   27.4  6.454244e-03 2.141655e-02
##   28.1  1.704022e-03 7.000467e-03
##   28.7  5.425739e-04 2.662738e-03
##   29.4  1.426235e-04 8.592894e-04
##   30    4.536407e-05 3.255965e-04
##   30.7  1.192024e-05 1.049064e-04
##   31.3  3.791117e-06 3.973184e-05
##   32    9.961548e-07 1.279902e-05
## 
##  Lower 95 Percent Confidence Limits
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3   High - 0     High - 1
##   0       0.9912555   0.9502673   0.7165756   0.8366981 0.99715478 9.996945e-01
##   0.653   0.9908775   0.9496746   0.7205324   0.8392318 0.99660966 9.996073e-01
##   1.31    0.9904778   0.9490632   0.7244073   0.8417018 0.99595527 9.994944e-01
##   1.96    0.9900620   0.9484425   0.7281329   0.8440663 0.99518309 9.993506e-01
##   2.61    0.9896245   0.9478050   0.7317479   0.8463507 0.99426285 9.991659e-01
##   3.27    0.9891567   0.9471393   0.7353016   0.8485870 0.99314769 9.989242e-01
##   3.92    0.9886710   0.9464645   0.7386825   0.8507058 0.99183702 9.986176e-01
##   4.57    0.9881589   0.9457692   0.7419412   0.8527396 0.99027491 9.982233e-01
##   5.22    0.9876185   0.9450521   0.7450728   0.8546865 0.98841323 9.977162e-01
##   5.88    0.9870386   0.9442998   0.7481175   0.8565718 0.98615752 9.970525e-01
##   6.53    0.9864345   0.9435333   0.7509773   0.8583356 0.98350730 9.962100e-01
##   7.18    0.9857952   0.9427393   0.7536930   0.8600040 0.98035072 9.951259e-01
##   7.84    0.9851071   0.9419028   0.7562959   0.8615967 0.97652946 9.937067e-01
##   8.49    0.9843881   0.9410463   0.7586989   0.8630612 0.97204557 9.919051e-01
##   9.14    0.9836246   0.9401546   0.7609340   0.8644176 0.96671406 9.895876e-01
##   9.8     0.9828000   0.9392098   0.7630212   0.8656786 0.96027466 9.865554e-01
##   10.4    0.9820036   0.9383131   0.7647492   0.8667176 0.95336712 9.830417e-01
##   11.1    0.9810131   0.9372170   0.7665473   0.8677926 0.94381724 9.777736e-01
##   11.8    0.9799506   0.9360613   0.7680929   0.8687096 0.93238126 9.708875e-01
##   12.4    0.9789770   0.9350180   0.7692005   0.8693606 0.92082640 9.633385e-01
##   13.1    0.9777606   0.9337320   0.7702175   0.8699494 0.90497463 9.520862e-01
##   13.7    0.9766426   0.9325641   0.7708325   0.8702956 0.88906659 9.398238e-01
##   14.4    0.9752412   0.9311152   0.7712210   0.8704969 0.86742606 9.216944e-01
##   15      0.9739487   0.9297899   0.7712440   0.8704788 0.84592741 9.021463e-01
##   15.7    0.9723229   0.9281332   0.7708687   0.8702106 0.81704140 8.736528e-01
##   16.3    0.9708179   0.9266054   0.7701635   0.8697450 0.78876238 8.434788e-01
##   17      0.9689176   0.9246780   0.7688371   0.8688921 0.75143047 8.005154e-01
##   17.6    0.9671517   0.9228830   0.7672141   0.8678612 0.71562740 7.563176e-01
##   18.3    0.9649131   0.9205943   0.7646734   0.8662585 0.66949610 6.956332e-01
##   18.9    0.9628244   0.9184380   0.7618631   0.8644918 0.62646612 6.358615e-01
##   19.6    0.9601659   0.9156535   0.7577311   0.8618979 0.57276614 5.579737e-01
##   20.2    0.9576754   0.9129939   0.7533448   0.8591434 0.52443219 4.857392e-01
##   20.9    0.9544923   0.9095072   0.7470733   0.8551976 0.46646243 3.979160e-01
##   21.6    0.9509737   0.9055171   0.7393365   0.8503119 0.40820365 3.102373e-01
##   22.2    0.9476572   0.9016048   0.7313237   0.8452262 0.35924119 2.388110e-01
##   22.9    0.9433935   0.8963268   0.7200574   0.8380248 0.30462141 1.642022e-01
##   23.5    0.9393606   0.8910428   0.7084547   0.8305407 0.26092759 1.109629e-01
##   24.2    0.9341585   0.8837522   0.6922050   0.8199334 0.21452202 6.395642e-02
##   24.8    0.9292231   0.8762793   0.6755256   0.8088816 0.17914515 3.672233e-02
##   25.5    0.9228400   0.8657091   0.6522578   0.7931644 0.14324933 1.760807e-02
##   26.1    0.9167706   0.8546003   0.6285129   0.7767364 0.11703505 8.815467e-03
##   26.8    0.9089073   0.8384952   0.5957030   0.7533319 0.09145063 3.725555e-03
##   27.4    0.9014216   0.8211872   0.5627122   0.7288949 0.07341441 1.721273e-03
##   28.1    0.8917183   0.7956362   0.5181564   0.6942884 0.05634446 6.798659e-04
##   28.7    0.8824824   0.7678786   0.4747835   0.6586120 0.04463258 3.013223e-04
##   29.4    0.8705217   0.7269205   0.4187943   0.6092226 0.03379979 1.148955e-04
##   30      0.8591572   0.6831003   0.3673832   0.5600263 0.02651346 4.980172e-05
##   30.7    0.8444797   0.6207078   0.3057136   0.4952240 0.01988446 1.862178e-05
##   31.3    0.8305846   0.5577061   0.2537627   0.4347508 0.01548829 7.968705e-06
##   32      0.8127241   0.4752420   0.1972006   0.3613996 0.01153507 2.945008e-06
##        group
## Day         High - 2     High - 3
##   0     1.000000e+00 9.997782e-01
##   0.653 1.000000e+00 9.997244e-01
##   1.31  1.000000e+00 9.996571e-01
##   1.96  1.000000e+00 9.995742e-01
##   2.61  1.000000e+00 9.994711e-01
##   3.27  1.000000e+00 9.993408e-01
##   3.92  9.999999e-01 9.991808e-01
##   4.57  9.999999e-01 9.989817e-01
##   5.22  9.999998e-01 9.987339e-01
##   5.88  9.999996e-01 9.984199e-01
##   6.53  9.999994e-01 9.980340e-01
##   7.18  9.999990e-01 9.975529e-01
##   7.84  9.999983e-01 9.969425e-01
##   8.49  9.999971e-01 9.961911e-01
##   9.14  9.999951e-01 9.952527e-01
##   9.8   9.999916e-01 9.940600e-01
##   10.4  9.999864e-01 9.927139e-01
##   11.1  9.999759e-01 9.907475e-01
##   11.8  9.999575e-01 9.882417e-01
##   12.4  9.999308e-01 9.855511e-01
##   13.1  9.998775e-01 9.816101e-01
##   13.7  9.998002e-01 9.773708e-01
##   14.4  9.996459e-01 9.711491e-01
##   15    9.994212e-01 9.644429e-01
##   15.7  9.989713e-01 9.545798e-01
##   16.3  9.983137e-01 9.439250e-01
##   17    9.969922e-01 9.282183e-01
##   17.6  9.950516e-01 9.112107e-01
##   18.3  9.911338e-01 8.860763e-01
##   18.9  9.853546e-01 8.587915e-01
##   19.6  9.736495e-01 8.183675e-01
##   20.2  9.563854e-01 7.743915e-01
##   20.9  9.216908e-01 7.091908e-01
##   21.6  8.607833e-01 6.251733e-01
##   22.2  7.761912e-01 5.351789e-01
##   22.9  6.268973e-01 4.081918e-01
##   23.5  4.568730e-01 2.863192e-01
##   24.2  2.459640e-01 1.516768e-01
##   24.8  1.076049e-01 6.851122e-02
##   25.5  2.816542e-02 1.995146e-02
##   26.1  6.914230e-03 5.580647e-03
##   26.8  1.119779e-03 1.061058e-03
##   27.4  2.150945e-04 2.319496e-04
##   28.1  2.957691e-05 3.659250e-05
##   28.7  5.233405e-06 7.208943e-06
##   29.4  6.779846e-07 1.048428e-06
##   30    1.160429e-07 1.968729e-07
##   30.7  1.464006e-08 2.751402e-08
##   31.3  2.465852e-09 5.039492e-09
##   32    3.068533e-10 6.891137e-10
## 
##  Upper 95 Percent Confidence Limits
##        group
## Day     Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3  High - 0   High - 1
##   0       0.9999862   1.0000000   1.0000000   1.0000000 0.9999765 1.00000000
##   0.653   0.9999839   1.0000000   1.0000000   1.0000000 0.9999685 1.00000000
##   1.31    0.9999811   1.0000000   1.0000000   1.0000000 0.9999577 1.00000000
##   1.96    0.9999779   1.0000000   1.0000000   1.0000000 0.9999433 1.00000000
##   2.61    0.9999742   1.0000000   1.0000000   1.0000000 0.9999241 1.00000000
##   3.27    0.9999699   1.0000000   1.0000000   1.0000000 0.9998980 1.00000000
##   3.92    0.9999648   1.0000000   1.0000000   1.0000000 0.9998635 1.00000000
##   4.57    0.9999590   1.0000000   1.0000000   1.0000000 0.9998174 1.00000000
##   5.22    0.9999521   1.0000000   1.0000000   1.0000000 0.9997559 1.00000000
##   5.88    0.9999441   1.0000000   1.0000000   1.0000000 0.9996724 1.00000000
##   6.53    0.9999349   1.0000000   1.0000000   1.0000000 0.9995623 1.00000000
##   7.18    0.9999242   1.0000000   1.0000000   1.0000000 0.9994157 1.00000000
##   7.84    0.9999115   1.0000000   1.0000000   1.0000000 0.9992168 1.00000000
##   8.49    0.9998971   1.0000000   1.0000000   1.0000000 0.9989556 1.00000000
##   9.14    0.9998805   1.0000000   1.0000000   1.0000000 0.9986084 1.00000000
##   9.8     0.9998608   1.0000000   1.0000000   1.0000000 0.9981389 1.00000000
##   10.4    0.9998403   1.0000000   1.0000000   1.0000000 0.9975780 1.00000000
##   11.1    0.9998127   0.9999999   1.0000000   1.0000000 0.9967102 0.99999999
##   11.8    0.9997805   0.9999999   1.0000000   1.0000000 0.9955379 0.99999997
##   12.4    0.9997488   0.9999999   1.0000000   1.0000000 0.9942131 0.99999993
##   13.1    0.9997062   0.9999998   0.9999999   0.9999999 0.9921764 0.99999981
##   13.7    0.9996643   0.9999996   0.9999999   0.9999999 0.9898867 0.99999956
##   14.4    0.9996083   0.9999994   0.9999998   0.9999999 0.9863886 0.99999879
##   15      0.9995533   0.9999991   0.9999997   0.9999998 0.9824843 0.99999713
##   15.7    0.9994801   0.9999985   0.9999995   0.9999996 0.9765718 0.99999216
##   16.3    0.9994085   0.9999977   0.9999992   0.9999995 0.9700414 0.99998153
##   17      0.9993135   0.9999963   0.9999986   0.9999991 0.9602770 0.99994999
##   17.6    0.9992210   0.9999944   0.9999979   0.9999986 0.9496521 0.99988304
##   18.3    0.9990988   0.9999909   0.9999965   0.9999978 0.9340511 0.99968684
##   18.9    0.9989804   0.9999863   0.9999946   0.9999966 0.9174305 0.99927675
##   19.6    0.9988247   0.9999778   0.9999911   0.9999946 0.8936350 0.99810093
##   20.2    0.9986746   0.9999667   0.9999864   0.9999919 0.8690091 0.99571307
##   20.9    0.9984784   0.9999465   0.9999780   0.9999870 0.8349231 0.98916357
##   21.6    0.9982572   0.9999146   0.9999645   0.9999794 0.7948669 0.97360319
##   22.2    0.9980460   0.9998729   0.9999469   0.9999696 0.7558840 0.94590269
##   22.9    0.9977724   0.9997991   0.9999157   0.9999526 0.7054766 0.88476971
##   23.5    0.9975130   0.9997044   0.9998759   0.9999310 0.6587542 0.79986403
##   24.2    0.9971792   0.9995397   0.9998072   0.9998944 0.6013319 0.66750510
##   24.8    0.9968649   0.9993329   0.9997218   0.9998493 0.5507629 0.54334440
##   25.5    0.9964634   0.9989832   0.9995794   0.9997750 0.4916777 0.41034263
##   26.1    0.9960882   0.9985581   0.9994086   0.9996867 0.4421170 0.31682661
##   26.8    0.9956125   0.9978680   0.9991356   0.9995466 0.3867976 0.23271938
##   27.4    0.9951711   0.9970677   0.9988235   0.9993874 0.3423158 0.17865853
##   28.1    0.9946157   0.9958412   0.9983504   0.9991471 0.2945147 0.13164000
##   28.7    0.9941042   0.9945086   0.9978393   0.9988882 0.2573598 0.10164795
##   29.4    0.9934651   0.9926123   0.9971108   0.9985194 0.2185890 0.07545664
##   30      0.9928805   0.9907089   0.9963726   0.9981452 0.1892137 0.05861793
##   30.7    0.9921551   0.9882143   0.9953875   0.9976443 0.1592146 0.04378889
##   31.3    0.9914957   0.9859022   0.9944520   0.9971666 0.1368990 0.03417847
##   32      0.9906823   0.9830877   0.9932803   0.9965651 0.1144526 0.02565550
##        group
## Day        High - 2  High - 3
##   0     1.000000000 1.0000000
##   0.653 1.000000000 1.0000000
##   1.31  1.000000000 1.0000000
##   1.96  1.000000000 1.0000000
##   2.61  1.000000000 1.0000000
##   3.27  1.000000000 1.0000000
##   3.92  1.000000000 1.0000000
##   4.57  1.000000000 1.0000000
##   5.22  1.000000000 1.0000000
##   5.88  1.000000000 1.0000000
##   6.53  1.000000000 1.0000000
##   7.18  1.000000000 1.0000000
##   7.84  1.000000000 1.0000000
##   8.49  1.000000000 1.0000000
##   9.14  1.000000000 1.0000000
##   9.8   1.000000000 1.0000000
##   10.4  1.000000000 1.0000000
##   11.1  1.000000000 1.0000000
##   11.8  1.000000000 1.0000000
##   12.4  1.000000000 1.0000000
##   13.1  1.000000000 1.0000000
##   13.7  1.000000000 1.0000000
##   14.4  1.000000000 1.0000000
##   15    1.000000000 1.0000000
##   15.7  1.000000000 1.0000000
##   16.3  1.000000000 1.0000000
##   17    1.000000000 1.0000000
##   17.6  1.000000000 1.0000000
##   18.3  0.999999998 1.0000000
##   18.9  0.999999987 1.0000000
##   19.6  0.999999898 0.9999999
##   20.2  0.999999404 0.9999995
##   20.9  0.999995367 0.9999963
##   21.6  0.999964753 0.9999755
##   22.2  0.999804571 0.9998821
##   22.9  0.998630262 0.9993197
##   23.5  0.993256843 0.9972467
##   24.2  0.963270221 0.9882810
##   24.8  0.877693124 0.9671076
##   25.5  0.673350712 0.9168243
##   26.1  0.464647518 0.8515175
##   26.8  0.271226227 0.7586888
##   27.4  0.163985742 0.6736794
##   28.1  0.089673379 0.5759389
##   28.7  0.053310204 0.4971782
##   29.4  0.029136942 0.4136558
##   30    0.017426496 0.3501580
##   30.7  0.009612636 0.2857524
##   31.3  0.005794909 0.2385446
##   32    0.003223455 0.1920658
eff.surv3 <- predictorEffect(c("Day"), survmod3a)
surv.effects3<-plot(eff.surv3,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "blue", "blue", "red", "red", "red", "red"), 
                              lty=c(1, 2, 4, 3, 1, 2, 4, 3)), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", 
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Fusion Partners")),
                                              cex=1, 
                                              cex.title=1)));surv.effects3

Display mean survivorship within multiple juvenile groups at the end of the experiment.

meanstresssurv2 <- ddply(multiple_final_stress, c("Fusion.Partners", "Treatment"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstresssurv2
##   Fusion.Partners Treatment  N       mean        sd         se max       lower
## 1               0   Ambient  9 0.55555556 0.5270463 0.17568209   1  0.02850928
## 2               0      High 26 0.07692308 0.2717465 0.05329387   1 -0.19482341
## 3               1   Ambient 16 0.81250000 0.4031129 0.10077822   1  0.40938711
## 4               1      High 15 0.06666667 0.2581989 0.06666667   1 -0.19153222
## 5               2   Ambient  9 0.77777778 0.4409586 0.14698618   1  0.33681923
## 6               2      High 16 0.00000000 0.0000000 0.00000000   0  0.00000000
## 7               3   Ambient 12 0.83333333 0.3892495 0.11236664   1  0.44408386
## 8               3      High  8 0.00000000 0.0000000 0.00000000   0  0.00000000
##       upper
## 1 1.0826018
## 2 0.3486696
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000

5. STRESS FUSION

Analysis

Analyze with a logistic binomial regression model.

Analyze fusion success over the course of the stress test.

#with full
fusemod2<-glmer(Fusion ~  Day * Treatment * Richness + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit, subset=c(Community=="Multiple")) 

summary(fusemod2)
Anova(fusemod2, type="II") 

plot(allEffects(fusemod2))

dispersion_glmer(fusemod2) #no overdispersion
plot(residuals(fusemod2)) + abline(h=0, lty=2)

meanstressfuse <- ddply(multiple_final_stress, c("Treatment", "Richness", "Timepoint"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressfuse
##   Treatment Richness Timepoint  N      mean        sd         se max
## 1   Ambient        1       End 19 0.7368421 0.4524139 0.10379087   1
## 2   Ambient        2       End 23 0.8260870 0.3875534 0.08081047   1
## 3   Ambient        4       End  4 1.0000000 0.0000000 0.00000000   1
## 4      High        1       End 25 0.6000000 0.5000000 0.10000000   1
## 5      High        2       End 18 0.5555556 0.5113100 0.12051692   1
## 6      High        3       End 14 0.4285714 0.5135526 0.13725270   1
## 7      High        4       End  8 1.0000000 0.0000000 0.00000000   1
##         lower    upper
## 1  0.28442818 1.189256
## 2  0.43853357 1.213640
## 3  1.00000000 1.000000
## 4  0.10000000 1.100000
## 5  0.04424556 1.066866
## 6 -0.08498116 0.942124
## 7  1.00000000 1.000000

Fusion is decreased at the end of the experiment.

6. STRESS GROWTH

Analyze with a poisson regression model.

Analyze growth during the course of the stress test.

Test for effect of temperature:

#with full dataset 
#polypmod2<-glmer(Polyps ~ Treatment + Community + Community:Treatment + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family = poisson, subset=c(Timepoint=="S2D11"))

#summary(polypmod2)
#Anova(polypmod2, type="II")
#qqPlot(residuals(polypmod2))

polypmod2a<-glmer(Polyps^2 ~ Day * Treatment * Community + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)

summary(polypmod2a)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Polyps^2 ~ Day * Treatment * Community + (1 | Parent.Site/Colony) +  
##     (1 | Treatment:Tank) + (1 | Slide/Sample)
##    Data: stress
## 
##      AIC      BIC   logLik deviance df.resid 
##   8676.7   8738.9  -4325.4   8650.7      870 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.9227 -0.8706  0.0316  0.8371  8.9788 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Sample:Slide       (Intercept) 0.198862 0.44594 
##  Slide              (Intercept) 0.004889 0.06992 
##  Colony:Parent.Site (Intercept) 0.191892 0.43805 
##  Treatment:Tank     (Intercept) 0.003544 0.05953 
##  Parent.Site        (Intercept) 0.011778 0.10853 
## Number of obs: 883, groups:  
## Sample:Slide, 156; Slide, 90; Colony:Parent.Site, 20; Treatment:Tank, 9; Parent.Site, 6
## 
## Fixed effects:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          3.498522   0.153886  22.735  < 2e-16 ***
## Day                                  0.010874   0.001113   9.773  < 2e-16 ***
## TreatmentHigh                       -0.158424   0.139109  -1.139   0.2548    
## CommunityMultiple                   -0.181614   0.124825  -1.455   0.1457    
## Day:TreatmentHigh                    0.014076   0.002115   6.654 2.85e-11 ***
## Day:CommunityMultiple                0.002843   0.001408   2.019   0.0435 *  
## TreatmentHigh:CommunityMultiple      0.265081   0.164704   1.609   0.1075    
## Day:TreatmentHigh:CommunityMultiple -0.031383   0.002530 -12.403  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Day    TrtmnH CmmntM Dy:TrH Dy:CmM TrH:CM
## Day         -0.108                                          
## TreatmntHgh -0.446  0.109                                   
## CmmntyMltpl -0.520  0.133  0.508                            
## Dy:TrtmntHg  0.049 -0.524 -0.159 -0.060                     
## Dy:CmmntyMl  0.085 -0.790 -0.086 -0.176  0.415              
## TrtmntHg:CM  0.343 -0.093 -0.770 -0.704  0.133  0.127       
## Dy:TrtmH:CM -0.041  0.438  0.133  0.089 -0.841 -0.555 -0.171
## convergence code: 0
## Model failed to converge with max|grad| = 0.00180889 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
Anova(polypmod2a, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps^2
##                            Chisq Df Pr(>Chisq)    
## Day                     302.2053  1  < 2.2e-16 ***
## Treatment                 0.8649  1    0.35238    
## Community                 4.2693  1    0.03881 *  
## Day:Treatment            48.8640  1  2.743e-12 ***
## Day:Community            34.2622  1  4.816e-09 ***
## Treatment:Community       0.2662  1    0.60586    
## Day:Treatment:Community 153.8404  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(residuals(polypmod2a))

## 246 982 
## 232 807

Plot difference in temperature between single and multiple groups:

stress$group<-paste(stress$Treatment, "-", stress$Community)

polypmod2b<-glmer(Polyps^2 ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)

eff.surv5 <- predictorEffect("Day", xlevels=4, polypmod2b)

polyp.effects2<-plot(eff.surv5,
                   lines=list(multiline=TRUE, 
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(1,2, 1, 2)), 
                   confint=list(style="bands", width=4), 
                   lwd=4,
                   #axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("(Total Polyp Expansion)"^2)), 
                   legend.position="right",
                   xlab=expression(bold("Time (Days)")), 
                   main="", 
                   lattice=list(key.args=list(space="top",
                                              border=FALSE, 
                                              title=expression(bold("Treatment - Number of Juveniles")),
                                              cex=1, 
                                              cex.title=1)));polyp.effects2

Present means of growth by number of juveniles. Sample size was too low due to mortality at the end of the experiment to calculate growth within multiple juvenile groups, so we pool these together for analysis here.

meanstressgrowth <- ddply(stress, c("Day", "Community", "Treatment"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressgrowth
##    Day  Community Treatment  N      mean       sd        se max     lower
## 1    0 Individual   Ambient 22  5.727273 2.333643 0.4975343  10  3.393630
## 2    0 Individual      High 23  5.000000 1.834022 0.3824200   9  3.165978
## 3    0   Multiple   Ambient 42  5.000000 1.847873 0.2851330   9  3.152127
## 4    0   Multiple      High 71  5.338028 1.764278 0.2093812   9  3.573750
## 5    3 Individual   Ambient 21  6.285714 2.052873 0.4479735  10  4.232842
## 6    3 Individual      High 24  5.833333 1.307725 0.2669383   9  4.525608
## 7    3   Multiple   Ambient 46  5.891304 1.636096 0.2412293   9  4.255208
## 8    3   Multiple      High 63  5.777778 1.754666 0.2210671  10  4.023112
## 9   11 Individual   Ambient 20  6.250000 2.022895 0.4523331  10  4.227105
## 10  11 Individual      High 19  6.105263 1.370107 0.3143241   9  4.735156
## 11  11   Multiple   Ambient 46  6.000000 1.646545 0.2427698  10  4.353455
## 12  11   Multiple      High 63  5.904762 1.672494 0.2107144  10  4.232268
## 13  17 Individual   Ambient 20  6.500000 2.139848 0.4784844  10  4.360152
## 14  17 Individual      High 18  6.166667 1.465285 0.3453709   9  4.701382
## 15  17   Multiple   Ambient 45  5.977778 1.815200 0.2705940  10  4.162578
## 16  17   Multiple      High 62  5.887097 1.747174 0.2218913  10  4.139923
## 17  23 Individual   Ambient 20  7.000000 2.294157 0.5129892  11  4.705843
## 18  23 Individual      High 16  6.250000 1.570563 0.3926406   9  4.679437
## 19  23   Multiple   Ambient 45  6.577778 1.888829 0.2815701  10  4.688948
## 20  23   Multiple      High 57  5.649123 1.894458 0.2509271  10  3.754664
## 21  26 Individual   Ambient 20  7.450000 2.480980 0.5547641  11  4.969020
## 22  26 Individual      High  4  8.250000 4.349329 2.1746647  12  3.900671
## 23  26   Multiple   Ambient 41  6.439024 2.013067 0.3143883  10  4.425957
## 24  26   Multiple      High 15  4.266667 1.869556 0.4827172   7  2.397111
## 25  32 Individual   Ambient 20  6.500000 2.781518 0.6219663  11  3.718482
## 26  32 Individual      High  2 12.000000 0.000000 0.0000000  12 12.000000
## 27  32   Multiple   Ambient 36  6.888889 2.213953 0.3689921  10  4.674936
## 28  32   Multiple      High  2  1.000000 0.000000 0.0000000   1  1.000000
##        upper
## 1   8.060915
## 2   6.834022
## 3   6.847873
## 4   7.102306
## 5   8.338587
## 6   7.141058
## 7   7.527401
## 8   7.532444
## 9   8.272895
## 10  7.475370
## 11  7.646545
## 12  7.577256
## 13  8.639848
## 14  7.631951
## 15  7.792977
## 16  7.634271
## 17  9.294157
## 18  7.820563
## 19  8.466607
## 20  7.543581
## 21  9.930980
## 22 12.599329
## 23  8.452091
## 24  6.136222
## 25  9.281518
## 26 12.000000
## 27  9.102842
## 28  1.000000
meanstressgrowth2 <- ddply(stress, c("Day", "Treatment"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressgrowth2
##    Day Treatment  N     mean       sd        se max    lower     upper
## 1    0   Ambient 64 5.250000 2.039296 0.2549121  10 3.210704  7.289296
## 2    0      High 94 5.255319 1.777616 0.1833471   9 3.477703  7.032935
## 3    3   Ambient 67 6.014925 1.770914 0.2163516  10 4.244011  7.785840
## 4    3      High 87 5.793103 1.636345 0.1754346  10 4.156758  7.429449
## 5   11   Ambient 66 6.075758 1.756838 0.2162518  10 4.318919  7.832596
## 6   11      High 82 5.951220 1.601715 0.1768799  10 4.349504  7.552935
## 7   17   Ambient 65 6.138462 1.919285 0.2380580  10 4.219177  8.057746
## 8   17      High 80 5.950000 1.683125 0.1881791  10 4.266875  7.633125
## 9   23   Ambient 65 6.707692 2.013417 0.2497336  11 4.694276  8.721109
## 10  23      High 73 5.780822 1.835200 0.2147940  10 3.945622  7.616022
## 11  26   Ambient 61 6.770492 2.209023 0.2828364  11 4.561469  8.979515
## 12  26      High 19 5.105263 2.941933 0.6749258  12 2.163330  8.047196
## 13  32   Ambient 56 6.750000 2.413974 0.3225809  11 4.336026  9.163974
## 14  32      High  4 6.500000 6.350853 3.1754265  12 0.149147 12.850853

Next, analyze the effect of temperature on growth (polyp expansion) at day 23 of 32 during the thermal stress test. As seen in table above, this is the timepoint with sufficient sample size, as later timepoints had high mortality.

First, for the effect of community type (individual vs multiple). No need to square transform as residuals are normal on scale of raw data.

polypmod3a<-glmer(Polyps ~ Treatment * Community + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=d23_stress, family = poisson)

summary(polypmod3a)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Polyps ~ Treatment * Community + (1 | Parent.Site/Colony) + (1 |  
##     Treatment:Tank) + (1 | Slide)
##    Data: d23_stress
## 
##      AIC      BIC   logLik deviance df.resid 
##    607.4    630.8   -295.7    591.4      131 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0539 -0.4361  0.0267  0.4877  1.6492 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 0.000e+00 0.000e+00
##  Colony:Parent.Site (Intercept) 8.515e-03 9.228e-02
##  Treatment:Tank     (Intercept) 1.351e-10 1.162e-05
##  Parent.Site        (Intercept) 0.000e+00 0.000e+00
## Number of obs: 139, groups:  
## Slide, 81; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      1.94217    0.08946  21.711   <2e-16 ***
## TreatmentHigh                   -0.11005    0.13286  -0.828    0.407    
## CommunityMultiple               -0.05117    0.10443  -0.490    0.624    
## TreatmentHigh:CommunityMultiple -0.05348    0.15657  -0.342    0.733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnH CmmntM
## TreatmntHgh -0.614              
## CmmntyMltpl -0.794  0.525       
## TrtmntHg:CM  0.523 -0.852 -0.665
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(polypmod3a, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                      Chisq Df Pr(>Chisq)  
## Treatment           4.5630  1    0.03267 *
## Community           0.9225  1    0.33683  
## Treatment:Community 0.1167  1    0.73266  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(residuals(polypmod3a))

## 380 265 
##  53  36

There was no difference in growth between individual and multiple juvenile groups, but an effect of temperature was significant. Next, analyze within multiple juvenile groups to look for effect of richness and fusion.

Generate a plot of this relationship.

d23_stress$group<-paste(d23_stress$Treatment, "-", d23_stress$Community)

polypmod3a2<-glmer(Polyps ~ Treatment + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=d23_stress, family = poisson)

eff.stress.polyps1 <- predictorEffect("Treatment", polypmod3a2)
polyps.stress.effects1<-plot(eff.stress.polyps1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), 
                   lwd=4,
                   axes=list(y=list(lim=c(5,8))),
                   type="response", 
                   ylab=expression(bold("Polyp Expansion")), 
                   legend.position="right",
                   xlab=expression(bold("Temperature")), 
                   main=""); polyps.stress.effects1

                   #lattice=list(key.args=list(space="right",
                                              #border=FALSE, 
                                              #title=expression(bold("Genotypic \nRichness")),
                                              #cex=1, 
                                              #cex.title=1)));surv.effects1
polypmod3b<-glmer(Polyps ~ Treatment * Richness * Fusion.Partners + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=d23_stress, family = poisson, subset=c(Community=="Multiple"))

summary(polypmod3b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Polyps ~ Treatment * Richness * Fusion.Partners + (1 | Parent.Site/Colony) +  
##     (1 | Treatment:Tank) + (1 | Slide)
##    Data: d23_stress
##  Subset: c(Community == "Multiple")
## 
##      AIC      BIC   logLik deviance df.resid 
##    458.4    490.0   -217.2    434.4       91 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.08863 -0.48696  0.01669  0.43919  1.42014 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 0.000e+00 0.000e+00
##  Colony:Parent.Site (Intercept) 1.293e-02 1.137e-01
##  Treatment:Tank     (Intercept) 0.000e+00 0.000e+00
##  Parent.Site        (Intercept) 6.070e-10 2.464e-05
## Number of obs: 103, groups:  
## Slide, 45; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                             1.61138    0.35341   4.560 5.13e-06 ***
## TreatmentHigh                           0.11566    0.41125   0.281    0.779    
## Richness                                0.13652    0.21716   0.629    0.530    
## Fusion.Partners                         0.15557    0.16543   0.940    0.347    
## TreatmentHigh:Richness                 -0.11739    0.23592  -0.498    0.619    
## TreatmentHigh:Fusion.Partners          -0.19732    0.21296  -0.927    0.354    
## Richness:Fusion.Partners               -0.06498    0.08519  -0.763    0.446    
## TreatmentHigh:Richness:Fusion.Partners  0.06960    0.09835   0.708    0.479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnH Rchnss Fsn.Pr TrtH:R TH:F.P Rc:F.P
## TreatmntHgh -0.840                                          
## Richness    -0.943  0.794                                   
## Fusn.Prtnrs -0.842  0.719  0.743                            
## TrtmntHgh:R  0.853 -0.932 -0.901 -0.675                     
## TrtmntH:F.P  0.656 -0.817 -0.582 -0.778  0.703              
## Rchnss:Fs.P  0.897 -0.759 -0.912 -0.911  0.824  0.711       
## TrtmH:R:F.P -0.776  0.871  0.790  0.787 -0.887 -0.904 -0.864
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(polypmod3b, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                                     Chisq Df Pr(>Chisq)  
## Treatment                          3.7823  1     0.0518 .
## Richness                           0.0517  1     0.8202  
## Fusion.Partners                    0.0047  1     0.9454  
## Treatment:Richness                 0.0793  1     0.7783  
## Treatment:Fusion.Partners          0.4499  1     0.5024  
## Richness:Fusion.Partners           0.0905  1     0.7635  
## Treatment:Richness:Fusion.Partners 0.5008  1     0.4792  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(residuals(polypmod3b))

## 380 862 
##  33  74

There is no effect of fusion or genotypic richness on growth.

Combine figures for Stress Period:

stress_figs<-plot_grid(surv.effects3b2, surv.effects3, polyps.stress.effects1, labels = c("A", "B", "C"), ncol=3, nrow=1, rel_heights= c(1,1,1), rel_widths = c(1,1,0.8), label_size = 20, label_y=1, align="h");stress_figs

ggsave(filename="Figures/stress_figs.png", plot=stress_figs, dpi=500, width=20, height=6, units="in")

7. TEMPERATURE DATA

Load in temperature data from Apex files. Data were recorded every 15 minutes and probes were calibrated to a digital thermometer weekly.

juvtemps<-read.csv("Data/ApexTemps.csv", header=T, na.strings="NA")
#convert Date and Time to date/time format
juvtemps$Date <- as.POSIXct(juvtemps$Date, format="%m/%d/%y %H:%M")

#convert to long format
juvtemps<-juvtemps %>% gather(Tank, Temperature, Tank17:Tank24)
juvtemps$Tank<-as.factor(juvtemps$Tank)

dailymax<-juvtemps[which(juvtemps$TimeDay == "930"),]
dailymax<-juvtemps[which(juvtemps$Phase == "Stress"),]

dailymax$Treatment<-ifelse(dailymax$Tank %in% c("Tank17"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(dailymax$Tank %in% c("Tank19"), "High",
                            ifelse(dailymax$Tank %in% c("Tank20"), "High",
                            ifelse(dailymax$Tank %in% c("Tank21"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank22"), "High",
                            ifelse(dailymax$Tank %in% c("Tank23"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank24"), "High",NA))))))))

maxtemps <- ddply(dailymax, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

juvtemps$Period<-ifelse(juvtemps$Phase %in% c("Stress"), "Stress",
                            ifelse(juvtemps$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtemps$Phase %in% c("Acclimation"), "Stress",
                            ifelse(juvtemps$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtemps$Phase %in% c("Ambient"), "Grow-Out",NA)))))

Plot timeseries of temperature data with tank temperature in colors.

#subset temperature measurements

juvtempsHIGH<-subset(juvtemps, Period==c("Stress"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHIGH$Treatment<-ifelse(juvtempsHIGH$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHIGH$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsHURR<-subset(juvtemps, Period==c("Hurricane"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHURR$Treatment<-ifelse(juvtempsHURR$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHURR$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsGROW<-subset(juvtemps, Period==c("Grow-Out"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsGROW$Treatment<-ifelse(juvtempsGROW$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsGROW$Tank %in% c("Tank19"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank20"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank22"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank24"), "Ambient",NA))))))))

juvtemps<-rbind(juvtempsHIGH, juvtempsGROW)
juvtemps<-rbind(juvtemps, juvtempsHURR)

TimeSeries1j<-ggplot(juvtemps, aes(x = Date, y = Temperature)) + 
  geom_line(aes(color = Treatment), size = 1) +
  scale_color_manual(values=c("blue", "red")) +
  ylab("Temperature (°C)")+
  xlab("Date")+
  ylim(26,32)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-07-31 10:00:00"))), linetype="solid", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-08-16 00:00:00"))), linetype="dashed", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-09-20 16:00:00"))), linetype="dotted", 
                color = "black", size=1)+
  theme_classic()+
  scale_x_datetime(limits = as.POSIXct(as.Date(c("2018-07-30 09:00:00", "2018-09-20 21:00:00")))) +
  theme(legend.position="none")+
  theme(text = element_text(size = 18, color="black"))+
  theme(axis.text = element_text(size=16, color="black"));TimeSeries1j

Plot by time of day by summarizing the mean temp at each time of day for each tank. Display the mean temperatures of each treatment during the Stress test phase (excluding the pause in treatment due to the hurricane system shutdown).

#subset juv temps for the stress period
juvtempsSTRESS<-juvtemps
juvtempsSTRESS$Period<-ifelse(juvtempsSTRESS$Phase %in% c("Stress"), "Stress",
                            ifelse(juvtempsSTRESS$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtempsSTRESS$Phase %in% c("Acclimation"), "Acclimation",
                            ifelse(juvtempsSTRESS$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtempsSTRESS$Phase %in% c("Ambient"), "Grow-Out",NA)))))

juvtempsSTRESS

juvtempsSTRESS<-juvtempsSTRESS[which(juvtempsSTRESS$Phase == "Stress"),]


QuarterlyJ <- ddply(juvtempsSTRESS, c("TimeDay", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

MeansJtemp <- ddply(juvtemps, c("Period", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )
MeansJtemp

MeansJtemp2 <- ddply(juvtemps, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   max = max(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )
MeansJtemp2

Generate final figure.

JuvTempsDay<-ggplot(QuarterlyJ, aes(x = TimeDay, y = mean, color=Treatment, fill=Treatment)) + 
  geom_line(color="black", size = 1) +
  #facet_wrap(~Period) +
  scale_x_continuous(breaks=c(0, 720, 1440), labels=c("0:00", "12:00", "24:00")) +
  scale_color_manual(values = c("blue", "red")) +
  scale_fill_manual(values = c("blue", "red")) +
  theme_classic() +
  geom_ribbon(aes(ymin=QuarterlyJ$lower, ymax=QuarterlyJ$upper), linetype=2, alpha=0.4, color=NA) +
  ylab("Temperature (°C)") +
  xlab("Time of Day") +
  ylim(26,32)+
  theme(text = element_text(size = 18, color="black"))+
  theme(panel.margin = unit(2, "lines"))+
  theme(axis.text = element_text(size=16, color="black")); JuvTempsDay

Combine figures for temperature.

temp_figs<-plot_grid(TimeSeries1j, JuvTempsDay, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(1,0.8), label_size = 20, label_y=1, align="h");temp_figs

ggsave(filename="Figures/temp_figs.png", plot=temp_figs, dpi=500, width=14, height=6, units="in")